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Abstract — A novel and computationally efficient algorithm for autonomous detection and localization of anomalies in 
hyperspectral imagery is presented. Anomaly refers to any object whose spectral radiance does not comport with that of its 
immediate neighborhood. It is assumed that the spatial extent of the anomaly is smaller than a sensor detector footprint, and 
that it is entirely confined to a single image pixel. The technique developed here is an unsupervised learning algorithm that 
examines each pixel in the context of its immediate neighborhood without any a priori knowledge about the spatial and 
spectral characteristics of the expected background or potential anomalies. The image representing each of the spectral 
bands of the hyperspectral image under consideration is independently converted to a two-dimensional binary anomaly map, 
which lends itself to straightforward parallelization of the computational process. The composite anomaly map is then 
obtained by adding the entire set of anomaly maps to which a threshold is applied and detection decisions are subsequently 
made. The results of the application of the algorithm to hyperspectral cubes obtained from the AVIRIS data and color RGB 
images are presented. It is shown that the algorithm provides a robust anomaly detection methodology with very-low 
computational overhead. This new algorithm has demonstrated computational efficiency of over three orders of magnitude 
better than the Boeing computationally -enhanced version of the N-FINDR. Unlike the N-FINDR, real-time application of the 
new anomalous source detection algorithm appears practicable. 
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I. Introduction 

The ability to detect potential threats in the early engagement phase of a mission improves the survivability of military 
assets and may provide enhanced protection for commercial aircraft. For example, the ability to detect an antiaircraft 
system while the physical distance between the aircraft and an unknown missile site is significant, offers greater 
opportunities for initiation of evasive action or deployment of defensive devices. The capability to scan a large area at high 
frame rate and detect objects that do not conform to their background is also of great value to many reconnaissance and 
surveillance systems. Long-range wide-area imagery of a suspected threat region can result in the concealed targets being 
confined to zones limited to less than a detector footprint. Algorithms for detection of targets smaller than the detector 
footprint cannot utilize spatial properties and must be entirely reliant on the spectral characteristics of potential targets . The 
objective of this research is to construct an algorithm to detect and locate subpixel anomalies in hyperspectral imagery 
(HSI). The algorithm must concurrently provide high probability of detection (PD) and low false -alarm rate (FAR) with 
the concomitant low-latency requirement of practical real-time systems. A computationally efficient and robust anomaly 
detection engine is a vital module in any hyperspectral target detection and tracking system where no a priori information 
about the spectral and spatial characteristics of the expected background and potential targets are available. Coupling such 
an algorithm with a fast frame -rate HSI data acquisition system, in principle, will make possible the simultaneous tracking 
of multiple targets each smaller than a detector footprint. 

The proliferation of HSI equipped sensor platforms, for both civilian and military applications, has led to the development 
of a plethora of analysis techniques for exploitation of the vast amounts of information contained within these images [1- 
12]. HSI comprises potentially hundreds of spatially co-registered images taken at narrow, and generally contiguous, 
spectral bands spanning visible, near infrared, mid- wave and long- wave infrared (IR), and/or millimeter wave of the 
electromagnetic spectrum. A HSI is a three-dimensional array consisting of one spectral and two spatial dimensions. Each 
spatially corresponding pixel of the data cube represents a portion of the spectrum of the emitted energy collected by the 
respective detector; this set of pixels along the spectral dimension is referred to as a spectral pixel array. The spectral 
radiance emanating from a pixel footprint containing an anomalous source is a blend of the anomalous source and a portion 
of the uncontaminated pixel spectral radiances. Theresulting signal from each detector in the spectral pixel array is 
determined by common radiometric methods [13-14]. The spectral signature of any pixel array is affected by the 
radiometric characteristics of materials located within the detector footprint, illumination, thermal self-emission, the 
intervening atmosphere, shadowing, and scattering. Spaceborne and airborne HSI, which are often used to provide detailed 
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information about the identity and relative abundance of earth surface materials, are widely used in environmental remote 
sensing and monitoring, commercial and agricultural land survey, and mineral resource exploration and assessment [3-8]. 
The utility of HSI imagery for detection, classification, identification and tracking of subpixel objects offers many practical 
applications in military, security, and commercial systems. 

Classification methods for HSI can be divided into two broad categories, viz., supervised and unsupervised. Supervised 
classification methods are based on the availability of a priori information regarding the spectral signatures of all surface 
materials of interest that may potentially be present in the image. Unsupervised classification procedures, however, start 
from a complete lack of a priori knowledge about the scene and anomalous sources that may potentially be present in the 
HSI and operate based upon a posteriori spectral information that is obtained from the input data cube. Extensive 
reflectance libraries of wide arrays of materials have been catalogued and numerous HSI algorithms for detection of objects 
of interest via exploitation of the known spectral signatures have been developed. Detection algorithms that utilize a priori 
knowledge about the spectral characteristics of expected target and background classes often rely on matched-filter (MF) 
theory [15]. 

Previous HSI anomaly detection algorithms are based on the theory of linear spectral unmixing, where the radiance of each 
pixel is assumed to be a weighted sum of the spectral signatures of a set of pure materials comprising the pixel area [16-21, 
31]. Over the years, many subtle variations of this theory have been formulated and an extensive assortment of 
computational techniques for the implementation of these formulations has been developed. Various versions of the linear 
mixing theory in their entirety share a set of common threads. For a typical approach, the minimum volume simplex in the 
spectral hyperspace that circumscribes the HSI data cube is constructed. The end-members, which are the vertices of the 
circumscribing simplex, constitute the spectral signatures associated with the set of pure materials comprising the entire 
image. Finear regression is applied on a pixel -by-pixel basis for spectral unmixing of the HSI in terms of the end- 
members. Two-dimensional abundance maps corresponding to each of the calculated end-members are computed, where 
the value of each pixel denotes the proportion of the pixel comprised of the respective end -member. End-members whose 
relative abundance throughout the HSI are minuscule in comparison to others are associated with anomalous materials, and 
pixels that include those end-members in their compositions are designated as anomalies. The anomaly detection 
algorithms based on spectral unmixing are plagued with extremely large computational overhead, making them unsuitable 
for real-time applications. 

The RX algorithm developed by I. S. Reed and X. Yu in 1990 and its many variations have gained wide acceptance and 
utility for HSI anomaly detection applications [22-28]. In general, the RX algorithm attempts to model the image 
background, which constitutes the preponderance of the materials comprising the HSI, as a set of multivariate Gaussian 
distributions in the spectral hyperspace. Each one of the HSI pixels is represented as a vector in the spectral hyperspace, 
and the entire set of pixels constitute one or more clusters with varied mass and shape in that space. Isolated points in the 
spectral hyperspace not contained within the clusters are associated with anomalous pixels. The metric utilized to quantify 
the degree of isolation of a pixel with respect to a particular cluster is the Mahalanobis distance. In most practical 
implementations of the RX algorithm, an annular region surrounding the pixel under consideration is used to model the 
local background. In order to obtain statistically reliable parameters for the Gaussian distributions representing local 
neighborhoods, the annular regions must include large numbers of pixels. This creates a problem when the set of pixels 
associated with the local area surrounding a potential anomaly have spectral signatures that are not tightly correlated to 
form a structurally cohesive cluster in the hyperspace. In such cases the algorithm generates a single multivariate Gaussian 
distribution for the local area which is not an accurate representation of the distribution. This results in computed 
covariance components that are far larger than the sum of respective covariance components of individual distributions, 
which in turn leads to significantly under-estimated point-cluster distance. For example, if the region surrounding an 
anomalous pixel comprises multiple disparate classes of spectral signatures, such is the case when the anomaly is located in 
the proximity of an edge, the anomaly will not be detected. 

The new HSI anomaly detection algorithm, presented in Section 2, is robust, conceptually intuitive, algorithmically simple, 
and computationally efficient. The performance of this algorithm is comparable to more exotic algorithms such as N- 
FINDR; a realistic comparison is included in Section 3.4. This algorithm can be implemented for real-time tracking of 
multiple moving subpixel targets in challenging environments using readily available processing platforms. In practical 
implementations, the subpixel anomaly detection step may be followed by deployment of an optical zooming mechanism to 
decrease the size of the detector footprint in the vicinity of the detected anomaly and allow further processing for target 
identification. Both the Boeing Enhanced N-FINDR and the new HSI anomaly detection algorithm fall under the 
unsupervised category. The significant difference between these methods and others is that the other methods attempt to 
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determine the abundance of specific materials in each pixel while the Boeing Enhanced N-FINDR and the new HSI 
anomaly detection algorithm attempt to locate pixels that are spectrally disparate from neighboring pixels while knowing 
nothing about the abundance or composition of any pixel. 

II. Subpixel Anomalous Source Detection Algorithm 

2.1 Concept 

The concept for the subpixel anomalous source detection (SASD) algorithm is based upon the idea that a detectable 
anomalous source located within a pixel footprint has a spectral radiance that differs from the spectral radiance of 
neighboring pixel footprints. 1 During this investigation, it was assumed that the anomalous source was totally contained 
within a pixel and that resultant radiance is computed by replacing a fraction of the original radiance with the same fraction 
of the anomalous source radiance as explained in Section 3.1. 2 This is referred to as pixel contamination. The process to 
be described converts an input three-dimensional (3D) data cube to a binary two-dimensional (2D) anomaly map as the 
output. The 3D data cube is generated from the sensor output and the data cube may represent RGB color, multispectral, 
hyperspectral, artificial color, or panchromatic grayscale images. The first two dimensions of the data cube represent 
spatial coordinates and the third denotes spectral image planes. Pixel values of unity in the output 2D anomaly map 
indicate input image coordinates containing anomalous sources, and zero -valued pixels are indicative of background 
regions. Generation of the anomaly map requires that each of the spectral image planes of the data cube under evaluation 
be first converted to three pseudo-images, viz., Laplacian, Edge, and Turbulence, which have the same spatial dimensions 
as the input image-plane. These are combined to form an incongruence map which is further transformed into the desired 
binary-valued anomaly map. Anomalies in the scene are searched for pixel -by-pixel while considering their surrounding 
neighborhood of pixels. In this investigation, the neighborhood comprises a 3x3 pixel array centered at pixel (m,n) as 
shown in Table 1. 


Table 1. Typical pixel located at (m,n) and its neighborhood. 

[m - 1 , n - 1 ) — (m-l,n + l) 

(m,n- 1 ) (m,n + 1 ) 

(m + l,n-l) (m + l,n) (m + l,n + l) 


2.2 Laplacian pseudo-image 

Applying the Laplacian filter (operator) to an input image is a well-known method for detection of abrupt signal variations. 
This filter is very effective for detection of localized anomalies such as signal bumps (abrupt signal increases) and dimples 
(abrupt signal decreases) in an input image plane of the 3D data cube. In regions of the input image where the signal is 
relatively constant or has a near-linear profile, the resultant Laplacian image has very small pixel values. This is a 
desirable property in that slowly varying background can be eliminated from the input image and allows concentration on 
potential regions of interest likely to contain anomalies. Portions of the image corresponding to shadow areas and regions 
approaching edges have, in general, relatively constant and near-bilinear intensities, respectively. In these image regions, 
the resultant Laplacian pseudo-image has near-zero pixel intensities. Each pixel of the Laplacian pseudo -image L is 

equal to the absolute value of the difference between the corresponding pixel in the source image -plane and the mean value 
of its eight neighbors, and is given by 


mnp 


n+\ m + 1 

-9D m „ n + X X D , , 
mnp , , t i mnp 

1 n=n - 1 m=m - 1 1 


0 < m < M -1, 0 < n< N - \, 1 <p<P 


( 1 ) 


where D and L denote, respectively, the source image -plane and its Laplacian pseudo-image, and M, A, P are the number of 
pixels along two spatial directions and the number of image planes in the 3D data cube, respectively. The value of each 
pixel in the Laplacian pseudo-image is proportional to the difference between the respective pixel of the input image and 


1 For simplicity, the pixel footprint of the detector imaged upon the ground by the sensor optics will generally be referred to as a pixel. 

2 In the remainder of this paper, the value associated with each pixel in the images comprising a hyperspectral data cube will be referred to as its signal 
rather than a radiometry quantity and each spectral image is considered as a grayscale image. Others factors such as detector linearity are ignored in this 
investigation. 
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the mean value of its close neighbors. However, the Laplacian operator has a serious disadvantage in that the resultant 
pseudo-image has large pixel values at image points corresponding to edges and naturally occurring segment boundaries. 

2.3 Edge pseudo -image 

In order to avoid problems associated with segment boundaries that plague the Laplacian pseudo -image, an Edge operator 
was developed during this investigation. Each pixel of the Edge pseudo -image is equal to the minimum of the absolute 
values of the difference between the corresponding pixel in the source image and its eight neighbors. The Edge pseudo - 
image is given by 

E mnp = \ D mnp ~ D m'n'p\’ ( 2 ) 

(m' ,n') 

m — 1 < in < m + 1 , n — 1 < li < n + \ 

where E denotes the Edge pseudo-image. This Edge filter tends to generate large output values for pixels in image areas 
approaching edges where signal varies linearly. However, it yields small output values for image pixels situated on the 
segment edges. This is due to the fact that edge pixels are connected to each other forming line segments and adjacent edge 
pixels have similar signal levels which result in lower values for output of the filter as stated in (2). 

2.4 Turbulence pseudo-image 

The signal at a particular pixel in the Turbulence pseudo-image is directly proportional to the degree of chaos associated 
with intensities of its eight neighboring pixels in the input image. The Turbulence pseudo-image is obtained by setting the 
value of each pixel equal to the unbiased estimator of the standard deviation of the eight neighbors of the respective pixel in 
the source image. 

1 n + 1 m + 1 / \2 

Tmnp = \ ^ ^ I ~ Emnp I '> ( m > n ) ^ ( m > w ) (3 a ) 

V 7 n =n - 1 m -m-1 \ 1 ’ 


and 
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(3b) 


where T represents the Turbulence pseudo-image andm is the neighborhood mean. It is noted that the unbiased variance 
is used in computation of the Turbulence pseudo-image. 

2.5 Incongruence and anomaly maps 

The method to declare the detection of potential anomalous sources requires that the anomaly occur in a user-specified 
minimum number of image planes/bands. Further, it does not matter which bands are designated for each anomalous 
source. Each source may have a different set of bands which may change over time due to temporal variation of view 
angles for example. To accomplish this, an incongruence map for each image plane is generated from which a binary 
anomaly map is produced that locates each pixel having been declared to contain an anomaly. For simplicity, the signal -to- 
noise ratio of each band was considered to be the same and was not included in the derivation of the incongruence map and 
anomaly map equations. 

The three non-negative pseudo images described by (1)— (3) are used to form P non-negative continuous -valued 
incongruence pseudo-images expressed by (4a), which is then threshold to form the binary incongruence map expressed by 
(4b). 
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where I and I denote, respectively, the continuous -valued incongruence pseudo-image and the binary- valued incongruence- 
map, and H is the user- specified incongruence threshold. Equations ( 1 )— (4) show that the process of obtaining the 
incongruence map corresponding to each one of the input source image planes is independent of the others and therefore 
can be carried out in parallel to maximize computational throughput. It is noted that when T mn p — » 0, the product 

hnnpEmnp ten ds to zero faster than T mn p thereby restricting I mnp from diverging. 


The non-negative integer-valued anomaly pseudo-image is then computed as the sum of the binary-valued incongruence 
maps corresponding to all the image planes/bands across the entire set of bands, which is then threshold to compute the 
binary anomaly map. Now, 

P 

A m n ~ X Imnp (^ a ) 

P = 1 


and 


1 ’ Ann - Q 
0; otherwise 


(5b) 


where A and A denote, respectively, the anomaly pseudo-image and anomaly map, and Q is the integer band-threshold 
parameter. The parameter Q is the aforementioned user-specified minimum-number of image bands necessary to declare 
pixel (m,n) as anomalous regardless of which mix of the P bands are used. 

Anomaly pseudo-image pixel values are directly proportional to values of the respective pixels in the Laplacian and 
Edgepseudo-images and are inversely proportional to the respective pixel value in the Turbulence pseudo -image. The 
value of the anomaly pseudo-image at a particular pixel is directly associated with the likelihood that the pixel under 
consideration is an anomaly. If a pixel of the anomaly pseudo-image has a relatively large value it implies that the pixel is 
not congruent with its neighborhood in the input image and is possibly an anomaly. Conversely, a small pixel value in the 
anomaly pseudo-image implies that the pixel in question is most likely a normal pixel in the input image. In the 
simulations presented in Section 3, it was found that for both of the ninety -band hyperspectral data cubes considered, 
setting incongruence and band thresholds to five and forty, respectively, resulted in satisfactory performance. Lowering the 
threshold values for H and Q , in general, leads to improved probability of detection at the expense of raising the false -alarm 
rate. For the large number of evaluations performed with RGB images, setting incongruence and band thresholds to, 
respectively, ten and two resulted in good performance. 

Table 2. Typical pixel coordinates comprising a new 2x2 pixel for each of the four similar images. 
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(m — l,n + 1)^ 
(m,n + 1) 


2.6 Expansion of methodology 

The anomalous source detection methodology presented makes the assumption that the anomalous source is confined 
within a single pixel. In practice, this likely will not be the situation since the subpixel -size anomalous source can straddle 
two, three, or four adjacent pixels. One relatively simple approach to resolve this possible difficulty is to artificially create 
new pixels comprising a 2 by 2 matrix of the original pixels, thereby forming a new image having one quarter of the 
original pixels. By shifting the reference pixel appropriately, three additional similar images are formed. The typical pixel 
coordinates comprising a new pixel for each of the four similar images are shown in Table 2. In a similar manner, such 
synthetic images can be formed having 3 by 3 matrix of the original pixels. The anomaly analysis would then be 
performed on each of these new images. 

An alternative approach is to enlarge the neighborhood surrounding the pixel and fit these neighboring pixels with, for 
example, a bicubic spline. A similar process to produce an anomaly map can be followed given appropriate modification to 
the equations. A further alternative approach is when multiple temporal frames of the scene are available and the sub-pixel 
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size anomaly is moving, it is believed possible to detect and track the anomaly even though the algorithm may not always 
be able to detect the anomalous source when it straddles multiple pixels in a given frame. 

III. Simulation Results 

In order to assess the efficacy of the subpixel anomalous source detection algorithm, a large number of evaluations were 
performed on actual images. Performance simulations were carried out on a wide selection of RGB images, a small subset 
of which is presented here, as well as two actual ninety-band hyperspectral data cubes. Numerous types of synthetic 
anomalies were transplanted into actual images in order to construct contaminated images constituting the simulated 
ground truth for each evaluation. The detection algorithm was applied to the contaminated evaluation images without any a 
priori knowledge about the image statistics and the locations and makeup of potential anomalies. The results pertaining to 
the coordinates of suspected anomalous pixels, as determined by the algorithm, were then compared to the ground truth for 
each evaluation and the detection and false-alarm rates were subsequently computed for each case. As will be observed 
from the following examples, the algorithm is capable of detecting subtle anomalies in complex imagery with high 
probability of detection and low false-alarm rate. For comparison purposes the performance of the Boeing Enhanced N- 
FINDR algorithm, a subpixel HSI anomalous -source detection method, is also presented in Section 3.4. 

3.1 HSI data preparation 

The two hyperspectral data cubes used to evaluate the algorithm were generated using AVIRIS data obtained from the 
NASA Jet Propulsion Laboratory [29]. One is a forest/rural scene from the “Low Altitude” data set and the other is a 
combined wilderness/urban scene from the “Jasper Ridge” data set, which will be henceforth referred to as Forest and 
Jasper-Ridge, respectively. Spatial dimensions of Forest and Jasper-Ridge images are 360x430 and 512x614 pixels, 
respectively. In order to keep the Boeing Enhanced N-FINDR processing within 2 GB of RAM when running in 
MATLAB, the 224 contiguous spectral bands with wavelengths from 400 to 2500 nanometers (nm) were down-sampled to 
90 channels using linear interpolation. The shortest wavelength channel was centered at 419.09 nm and the longest 
wavelength channel was centered at 2448.26 nm, with a bandwidth of 22.80 nm for each downsampled channel. In order 
to synthesize anomalies, a camouflage paint spectrum, henceforth referred to as “Camopaint,” was also generated from 
AVIRIS data of an olive-green paint spectrum that was likewise obtained from the NASA Jet Propulsion Laboratory. It 
was downsampled by linear interpolation to have the same spectral channel bandwidths as the two aforementioned 
hyperspectral data cubes. The Camopaint was selected because its spectrum is generally similar to the spectral content of 
the images and therefore should be more stressing than if something like aluminum or blue paint was selected as the 
anomaly. 



50 100 150 200 250 300 350 400 


Figure 1. Grayscale image representation of the ninety-band Forest data cube. 


Figure 1 shows the grayscale image representation of the Forest evaluation image, which was obtained by the pixel -wise 
summation of the entire set of ninety image planes of the respective down-sampled HSI data cube. Plots in Fig. 2 show 
spectral signatures of eight pixels that illustrate the spectral diversity in the Forest data cube. In these plots, the abscissa 
and ordinate axes represent, respectively, band number and relative signal. 


Page I 57 





International Journal of Engineering Research & Science (IJOER) 


ISSN: [2395-6992] 


[Vol-3, Issue-2, February- 2017] 




Figure 2. Spectral signatures of eight representative pixels chosen from eight 

DIFFERENT REGIONS OF FOREST DATA CUBE. 


Since the center wavelength and full width of the n th channel of any pixel of Forest and Jasper-Ridge data cubes and those 
of the Camopaint spectrum are identical, a weighted spectrum of any combination of image pixels and Camopaint is easily 
obtained on a channel-by-channel basis. The spectral signature of Camopaint, which is used as the source of pixel 
contaminations in the following examples, is shown in Fig. 3. 


In the subsequent examples, the input HSI evaluation images, in each case, were obtained from data cubes Forest and 
Jasper-Ridge whose elements are /(ra,rc,p), where /represents the unadulterated, or original, HSI data cube with m and n 
denoting spatial coordinates, and p is the band number. Anomalies are synthesized by mixing one or more pixels of the 
original HSI with Camopaint represented by spectral signature c (p) as shown in Fig. 3, where P denotes the total number of 
bands. The contaminated image is obtained as follows. 


f(m,n,p) 


|(1 -R)f(m,n,p) + (a)(R)c(p); ( m,n ) = ( m',n' ) 


( 6 ) 


where 


a= £ f{m\n\p) / Z c(p) (7) 

P = 1 / P= l 


and f(m,n,p) represents the contaminated hyperspectral image synthesized to be used as the input evaluation image with 
(m\ n') denoting the spatial coordinates of the contaminated pixel. The contamination factor R represents the fraction of 
Camopaint content in the contaminated pixel, where R = 1 indicates total replacement of the original signal with Camopaint 
and R = 0 means the original content of the pixel remains intact. The spectral profile of the contaminated pixel ( m' , n') in 
the synthesized contaminated data cube is obtained by mixing, in accordance with (6), the original content of this pixel with 
Camopaint. For the purpose of this investigation, the stressing case is considered to occur when the contaminated pixel is 
constructed such that the total or spectrally-integrated signal is equal to the original uncontaminated pixel spectrally- 
integrated signal. Although there are alternative ways to contaminate the pixel, the procedure used equates the spectral 
signal of the original pixel / with that of a fully contaminated pixel as expressed by (6) where a represents the ratio of the 
spectrally-integrated signal of the uncontaminated pixel to the spectrally -integrated value of the Camopaint. Clearly, this 
construction of a contaminated pixel maintains the total integrated signal of the modified pixel at a constant value 
regardless of the value of R where 0 < R < 1 . If now (7) is generalized so that spectral signature of the contaminant can be 
pixel dependent, then (8) expresses the constant composite signal. 

So, 

P p 

H f(m,n,p) = I.hl-R m , n )f(m,n,p) + a mn R mn c mn (p)] (8) 

P= 1 p= 1 L J 

where R mn , cc mn , and c mn (p ) refer to the values for pixel ( m,n ). In actual situations, the contamination could be 

significantly greater or smaller than that of the original pixel; however, the aforementioned contaminated pixel construction 
was used to explore the marginal detection capability of the subpixel anomalous source detection algorithm. 
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Figure 3. Spectral signature c(p) of Camopaint where p is the band number. 




(a) (b) 

Figure 4. (a) Grayscale image representation of the ninety-band Jasper-Ridge data 
cube, (b) Spectral signatures of five diverse pixels of Jasper-Ridge chosen from five 

DIFFERENT REGIONS, WHERE PIXEL COORDINATES ARE SHOWN IN THE LEGEND 


Figure 4(a) shows the grayscale image representation of the Jasper-Ridge HSI, which was obtained by the pixel -wise 
summation of all the ninety image planes of the respective HSI data cube. Plots of the spectral signatures of five diverse 
pixels from the Jasper-Ridge HSI are shown in Fig. 4(b). 


3.2 Performance results for hyperspectral evaluation images 

In the following examples, the hyperspectral evaluation images are constructed by creating synthetic signal functions and 
transplanting them in the pixels of Forest and Jasper-Ridge data cubes. Pixel coordinates are chosen randomly and their 
original contents are replaced with synthetic transplants. For each evaluation, the contamination factor R is set to a fixed 
value and each randomly chosen pixel is replaced with a transplant that is synthesized as thelinear mixture of the original 
content of that pixel and Camopaint in accordance with (6). The spatial coordinates of the transplanted pixels constitute the 
ground truth for each evaluation image. For each set of user-prescribed parametric values, Monte Carlo simulations were 
performed by choosing one-hundred pixel positions at random and replacing them with synthetic signal functions as 
previously described in order to construct the input HSI evaluation image. The algorithm is then applied to the HSI 
evaluation image and the anomaly-map is computed in accordance with the procedure of Section 2. The spatial coordinates 
of anomalies determined from examination of the output anomaly-map are then compared to the ground truth, and 
detection and false-alarm rates are computed. Detection rate is the percentage of anomalous pixels in the HSI evaluation 
image that are correctly identified, and false-alarm rate is the proportion of the number of pixels that are erroneously 
misclassified as anomalous. Each evaluation was repeated ten times and the detection and false-alarm rates for each case 
were computed and subsequently averaged across all trials. The entire process is then repeated using a different 
contamination factor R. A quite severe criterion was used for declaring a potential anomalous source a false alarm in that if 
the actual anomalous source was not designated in the correct pixel location (ra', n'), then it was declared a false alarm. 
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Figure 5. Anomaly detection rate as 

FUNCTION OF PIXEL CONTAMINATION FACTOR R 

for Forest. Band-threshold is set at Q = 

30, AND THE FUNCTION IS PLOTTED FOR THREE 
DIFFERENT VALUES OF INCONGRUENCE- 
THRESHOLD H = 3, 4, and 5. Average false- 
ALARM RATES PER ONE-MILLION INPUT PIXELS 
ARE 37.88, 6.46, AND 0, RESPECTIVELY 
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Figure 6. Anomaly detection rate as 

FUNCTION OF PIXEL CONTAMINATION FACTOR R 

for Forest. Incongruence-threshold is set 
at H = 5. Each plot corresponds to a fixed 
VALUE BAND-THRESHOLD Q = 10, 20, 30, AND 40. 
FALSE-ALARM RATES ARE ZERO FOR ALL FOUR 
CASES. 


Simulation results using HSI evaluation images based on Forest are presented in Figs. 5 and 6 which provide plots of 
detection rate as function of pixel contamination factor R for fixed values of band and incongruence thresholds. In Fig. 5, 
the value of band-threshold was set at Q- 30. The contamination factor R was varied from 0.1 to 1 for three different 
settings of the incongruence-threshold, viz., H- 3,4, and 5. It is evident that as the incongruence-threshold H increases, the 
average false-alarm rate decreases. In Fig. 6, the value of incongruence-threshold was set at H = 5 and the contamination 
factor Rwas varied from 0.1 to 1 for four different settings of the band-threshold, viz., Q = 10, 20, 30, and 40. As expected, 
increasing the contamination factor resulted in better detection performance with virtually no effect on the false-alarm rate. 
Raising the values of incongruence and band thresholds leads to lower detection and false-alarm rates. The plots in Figs. 5 
and 6 show that the algorithm performs remarkably well since it detects the anomalies in the HSI evaluation image with 
high probability of detection and very low false-alarm rate. For example, having a pixel contamination factor of R = 1, 
with incongruence and band threshold values of H = 5 and Q = 30, resulted in detection of all one-thousand transplanted 
anomalies (ten trials of one-hundred randomly placed transplanted anomalies in each case) with zero false-alarm. Even for 
pixel contamination factor as low as R = 0.5, close to ninety -percent of anomalies were detected with zero false-alarm. 
This is an impressive feat considering that each transplanted pixel preserves fifty -percent of its original content and yet the 
algorithm correctly isolates ninety -percent of the one-thousand contaminated pixels without misclassifying a single 
unadulterated Forest pixel. Figure 5 also shows that one can obtain higher detection rate at the expense of higher false- 
alarm rate. For example, lowering the value of incongruence -threshold from five to four resulted in higher probability of 
detection at the expense of raising the number of false alarms from zero to six pixels per million pixels in the HSI 
evaluation image. 

Figures 7-9 represent the results of simulations performed on the Jasper-Ridge HSI. As before, each evaluation involves 
the random selection of one-hundred pixels and replacing the content of each pixel with the synthetic signal function 
obtained by linearly mixing the original content and Camopaint. Each simulation was repeated ten times and average 
detection and false-alarm rates were computed. In the example of Fig. 7, the incongruence-threshold was set at H = 10, and 
the detection rate was computed as function of contamination factor R for four values of band-threshold, viz., Q= 10, 20, 
30, and 40. As expected, higher contamination factor and lower band-threshold both result in improved detection. Figure 8 
shows parametric plots of the number of misclassified pixels per one-million input-pixels as function of the incongruence 
threshold for a fixed band-threshold Q = 40. Plots are shown for four values of the contamination factor R = 0.10, 0.25, 
0.50, and 1. It is observed that raising the value of incongruence -threshold leads to lower false-alarm rates. The value of 
the contamination factor has virtually no effect on false-alarm rate and the variations in Fig. 8 are due to the statistical 
nature of the problem. Plots in Figs. 7 and 8 show that for a contamination factor of R = 1 with parametric settings of Q = 
40 and H - 10 result in detection of all the anomalous pixels with virtually a zero false-alarm rate. Even for a 
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contamination factor of R = 0.75, it is possible to achieve eighty-percent detection rate with zero false-alarm rate. Figure 9 
shows the parametric plots of the detection rate as function of incongruence threshold for a fixed band-threshold. Figures 8 
and 9 show that setting Q - 40 and H - 5 results in detection and false-alarm rates of one-hundred percent and zero, 
respectively, for a contamination factor of R = 1. Even for a contamination factor as low as R = 0.5, a detection rate of 
better than ninety -percent concurrent with zero false-alarm rate was achieved. 
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Figure 7. Anomaly detection rate as 

FUNCTION OF PIXEL CONTAMINATION FACTOR R 

for Jasper-Ridge. Incongruence- 
threshold IS SET AT H = 10, AND THE 
FUNCTION IS PLOTTED FOR FOUR VALUES OF 
BAND-THRESHOLD Q = 10, 20, 30, AND 40. 
False-alarm rates are zero for all four 

CASES. 


Figure 8. Number of falsely labeled 

PIXELS PER ONE-MILLION INPUT PIXELS IS 
PLOTTED AS FUNCTION OF INCONGRUENCE- 
THRESHOLD H. Band-threshold is fixed at Q 
= 40 AND FALSE-ALARM PERFORMANCE IS 
PLOTTED FOR CONTAMINATION FACTOR (CF) OF 
R = 0.25, 0.5, 0.75, AND 1. 


forty bands 



threshold 


Figure 9. Anomaly detection rate versus incongruence-threshold H with band 

THRESHOLD SET AT Q = 40. PIXEL CONTAMINATION FACTOR R IS THE CONTROL PARAMETER. 

3.3 Performance results for color RGB images 

A quantitative assessment of the SASD algorithm was performance as it applies to RGB images. Fig. 10(a) shows the 
RGB evaluation image in its original form, and Fig. 10(b) shows the RGB image used for sources of pixel contamination. 
Ten contaminant RGB vectors were synthesized, each by computing the mean value of the RGB vectors of all the pixels in 
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one of ten rectangular regions similar to the three marked in red in Fig. 10(b). This was done in order to emulate the effect 
of increased range on the synthesized tank-RGB -vectors that are used as contaminants. One-hundred pixels of the image in 
Fig. 10(a) were selected randomly and the RGB vector of each selected pixel was replaced with one of the ten synthesized 
tank-RGB -vectors, chosen randomly, in order to construct the evaluation image. The anomaly detection algorithm was 
then applied to the evaluation image. Detection and false alarm rates were computed by comparing the anomaly 
coordinates indicated by the output anomaly map with the actual positions of the transplanted pixels. The simulation was 
repeated one-hundred times and the results were averaged. It is noted that pixel contamination factor for this example was 
R = 1. Plots in Figs. 11 and 12 present the detection and false-alarm rates as functions of the incongruence threshold H 
with band threshold Q as the parameter. With incongruence and band thresholds set to 4 and 2, respectively, ninety -percent 
detection rate and zero false-alarm were observed. 

3.4 Enhanced N-FINDR Comparison 

The Boeing Company developed an Enhanced N-FINDR algorithm which follows the general framework of hyperspectral 
unmixing techniques such as N-FINDR [18, 30]. Both of these techniques typically determine the basic material 
components, or “endmembers,” of a hyperspectral scene without regard to the spatial correlation of the pixels. Instead, the 
scene pixels are treated as a statistical ensemble of spectra where each spectrum (pixel) is assumed to be a weighted sum of 
the scene endmember spectra. Statistical techniques are used to determine an optimal set of endmember spectra 
corresponding to the ensemble. Of course, a truly anomalous pixel would be expected to have endmember content 
substantially different from the overwhelming majority of proximate pixels in the scene. 



(a) (b) 

Figure 10. (a) RGB image to be corrupted by randomly selecting some of its pixels and 

REPLACING THEM WITH THE CONTAMINANT RGB VECTORS OBTAINED FROM THE IMAGE IN FlG. 10(B). (B) 

Contaminant RGB vectors are obtained by averaging the RGB vectors of ten rectangular 

PATCHES SIMILAR TO THE THREE MARKED REGIONS OF THE RGB IMAGE SHOWN. 
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Figure 11. Anomaly detection rate versus 

INCONGRUENCE THRESHOLD H WITH BAND THRESHOLD 
Q USED AS THE PARAMETER. 



Figure 12. False-alarm rate versus 

INCONGRUENCE THRESHOLD H WITH BAND 
THRESHOLD Q USED AS PARAMETER. 
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Figure 13. Enhanced N-FINDR results for 
Forest 


Enhanced N-FINDR with residual threshold = 0.8% 



Figure 14. Enhanced N-FINDR results for 
Jasper-Ridge 


Subpixel anomaly detection evaluations were performed on the Jasper-Ridge and Forest data cubes contaminated with 
Camopaint spectra using both the Boeing Enhanced N-FINDR anomaly detection engine and the SASD algorithm. Figures 
13 and 14 present the detection rate versus contamination factor for Forest and Jasper-Ridge, respectively. To compare 
these cases having similar detection performance for small contamination factor, the Boeing algorithm thresholds were set 
to give detection results similar to the subpixel anomalous source detection algorithm with the contamination factor R = 
0.2. The Boeing algorithm is conceptually similar to the original N-FINDR, but is far more computationally efficient. 
Accordingly, this represents a hyperspectral unmixing approach with minimal computational overheard. Nevertheless, the 
unmixing operation required approximately 630 MB of RAM and 8. 1x10 10 multiplications and additions for the Jasper- 
Ridge scene and approximately 300 MB of RAM and 2.6 x 10 10 multiplications and additions for the Forest scene. On the 
other hand, the SASD algorithm required only 1.54 x 10 7 additions and multiplications for the Jasper-Ridge scene and 7.6 x 
10 6 additions and multiplications for Forest scene. Although the Boeing Enhanced N-FINDR yields better detection rates 
than the SASD algorithm for contamination factors greater than R = 0.3, it achieves these results at the expense of 
dramatically greater computational overhead; however, the false-alarm rates are greater for the Boeing Enhanced N-FINDR 
than for the SASD algorithm when R> 0.3. For applications that cannot tolerate the additional overhead, the SASD method 
offers an alternative with an impressive balance between detection and false alarm rates. In these particular evaluations, the 
SASD algorithm required about 0.02% of the computations required by the Boeing Enhanced N-FINDR. For applications 
demanding higher detection performance, perhaps using a hybrid of the SASD technique and the Boeing hyperspectral 
unmixing technique can provide satisfactory detection performance with significantly reduced computational overhead. 

IV. Conclusions 

A novel method has been presented for detection and localization of subpixel anomalies in hyperspectral imagery. The 
conceptually intuitive HSI anomalous source detection algorithm, SASD, has been shown to be robust and remarkably 
computationally efficient. This anomaly detection procedure does not rely on any a priori information about the spatial and 
spectral characteristics of the expected background or of potential anomalies. The subpixel anomalous source detection 
algorithm converts the input HSI data cube into a 2D binary anomaly map, using two user-prescribed thresholds, in a 
computationally straightforward approach. Each pixel is examined in the context of its immediate neighborhood in each of 
the HSI image planes and an incongruence score is assigned to it which results in an incongruence pseudo-image for each 
spectral band. The incongruence score of a pixel has an inverse relationship to the degree of its conformity to its 
neighborhood. The user-prescribed incongruence threshold is used to convert each of the band-specific incongruence 
pseudo-images into a respective binary incongruence -map. These maps are subsequently combined to compute the binary 
anomaly map using the user-prescribed band-threshold. It has been shown through many simulations, using actual HSI and 
RGB images, that the subpixel anomalous source detection algorithm (SASD) is capable of detecting subtle anomalies with 
high probability of detection and low false-alarm rates. Comparisons with the Boeing Enhanced N-FINDR algorithm has 
confirmed that the computational processing power required for the implementation of the described subpixel anomalous 
source detection algorithm is several orders of magnitude less than for the Enhanced N -FINDR. The SASD algorithm can 
be implemented for processing high spatial-spectral resolution HSI with low latency requirements using commercial off- 
the-shelf (COTS) processing platforms. The low computational overhead afforded by this algorithm makes it a good 
candidate for real-time processing of high frame -rate HSI video streams in both commercial and military applications. 
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Future research is planned to (1) further understand the relationship of hyperspectral image data cubes' complexity and 
unresolved and subpixel-scale object anomalies detect ability from their immediate neighborhood, (2) investigate spectral 
bands selection of image data cubes to improve performance and reduce computational resources, (3) study the efficacy of 
algorithms if synthetic resolution-enhancement is applied to image data cube, (4) experimentally validate the algorithms, 
and (5) extending the algorithm to include imagery signal-to-noise information. It is anticipated that providing technology 
to detect and locate unknown unresolved and subpixel -scale anomalous objects in unknown scenes can positively impact 
industrial and military security needs by affording the ability to detect potential threats in the early engagement phase of a 
mission to improve (i) survivability of military assets and enhance protection for commercial aircraft, (ii) locating buried 
mines and IED detection, (iii) monitoring of pollution in smokestack plumes, (iv) chemical and biological agent detection, 
and (v) detection of small floating objects to aid in search and rescue operations to name a few potential applications. 
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